Title: “Modeling Coral Invertebrate Communities Using 3D photogrammetry”

author: “Journ Galvan” date: “4/23/2020” output: html_document —

# Load libraries
library(here)
library(dplyr)
library(tidyverse)
library(vegan)
library(reshape)
library(ggplot2)
library(ggpubr)
library(zoo)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(psych)
library(knitr)
library(faraway)
library(car)
library(MASS)
library(gridExtra)
library(grid)
library(AER)

#Load data 
branch_width <- read.csv(here("photogrammetry_data/branchwidth_data.csv")) 
coral_pg <- read.csv(here("photogrammetry_data/galvan_journ_datasheet_v2.csv")) 
coral_field <- read.csv(here("photogrammetry_data/field_experiment_colony_measurements_moorea_summer2019.csv"))
cafi_surveys <- read.csv(here("cafi_data/revised_cafi_data_moorea_summer2019_11_27.csv"))
cafi_field <- read.csv(here("cafi_data/prelim_cafi_counts_moorea_summer2019.csv"))
updated_cafi <- read.csv(here("cafi_data/cafi_data_w_taxonomy_summer2019_2020_5_21.csv"))

Incorporate cleaned data taken in the field of coral associated fish and invertebrates

Data was filtered to only include invertebrates identified in coral colonies. These were identified down to the species level. Species richness, abundance and Shannon weiner diversity index was calculated.

#Incorporate cleaned CAFI data
updated_cafi2 <- updated_cafi %>% filter(str_detect(coral_id, "^FE")) %>%
  dplyr::select(master_sort, coral_id, code, type, search_term, lowest_level, phylum, genus, species, general_notes) %>%
  subset(phylum!="Chordata") #filter out vertebrates

#Calculate CAFI richness, abundance, and diversity (shannon weiner) for each coral 
cafi_summarized2 <- group_by(updated_cafi2, coral_id) %>%
  summarise(num_cafi = n(), cafi_richness = length(unique(code)), cafi_present = paste(sort(unique(code)), collapse = ";"))
cafi_summarized2$sw <- updated_cafi2 %>% 
  count(code, coral_id = coral_id) %>% 
  spread(code,n) %>% 
  mutate_all(list(~tidyr::replace_na(.,0))) %>% 
  dplyr::select(-coral_id) %>% 
  diversity(index = "shannon")

Incorporate coral data

Cleaned data of coral colonies taken from the field and using photogrammetry techniques were combined into one dataframe. Many more and accurate morphological measurements could be acquired using photogrammetry than with traditional measurements taken *in situ. Common measurements of volume, surface area and available space known as interstitial space were calculated of each coral colonly.

Available space was traditionally measured using 2D measurements of five randomly selected branch distances and averaging them. In this experiment, we took methods used by Neil E. Doszpot et al. to calculate a 3D measurement of available space using convex hull geometry and the software estimated coral volume.

Because we did not take random branch distance measurements in the field, we used the software models to randomly select branch distances to be averaged.

Three other measurements were calculated from the coral models: convexity and sphericity which capture volume compactness and packing which captures how much of an objects surface area is situated internally versus externally.

Wide or tight branching corals were initially qualitatively determined. However, classifying corals with convexity values \(>= 0.5\) as “tight” and \(< 0.5\) as “wide” offers another method of classifying coral morphology. For the rest of these tests, we will classify “wide” and “tight” branching corals using this method.

# Clean field data
coral_field2 <- coral_field %>%
  dplyr::rename(branch = branch_width)

# Clean photogrammetry morphometric data
coral_pg$volume_pg <- as.numeric(as.character(coral_pg$volume_pg)) # Convert factor to numeric

# Clean photogrammetry branch width data and take averages of branch distance for each coral
branch_width$branch_distance_mm <- as.numeric(as.character(branch_width$branch_distance_mm)) # Convert factor to numeric
branch_w_summarized <- group_by(branch_width, coral_id) %>%
  summarise(avg_w_cm = sum(branch_distance_mm/10)/n(), #take average branch distance and convert mm to cm
            measurements = length(unique(replicate_measurement)), 
            locations = paste(sort(unique(location)), collapse = ";"))

#Create dataframe called coral_dim 
coral_dim <- merge(coral_pg, branch_w_summarized, by = "coral_id") %>%
 merge(coral_field2, by = "coral_id") %>%
  mutate(volume_pg=volume_pg*10^6,
         max_hull_volume=max_hull_volume*10^6, #convert m^3 to cm^3
         max_hull_surface_area=max_hull_surface_area*10^4,
         surface_area=surface_area*10^4,
         height_pg=height_pg*100,
         length_pg=length_pg*100,
         width_pg=width_pg*100) %>%   #convert m^2 to cm^2
  dplyr::rename(max_hull_SA = max_hull_surface_area,
                SA = surface_area)
          
coral_dim$interstitial_space <- coral_dim$max_hull_volume - coral_dim$volume_pg #calculate available space by subtracting software estimated volume from convex hull volume

coral_dim$SAV <- coral_dim$SA / coral_dim$volume_pg #calculate surface area to volume relationship

coral_dim$convexity <- coral_dim$volume_pg / coral_dim$max_hull_volume #calculate proportion occupied, high ratios indicate less free space between branches

coral_dim$packing <- coral_dim$SA / coral_dim$max_hull_SA #how much of an objects surface area is situated internally

coral_dim$sphericity <- ((3.14^(1/3))*((6*coral_dim$volume_pg)^(2/3)))/coral_dim$SA #calculate sphericity or how close the object is to a sphere

coral_dim$diff_est <- coral_dim$volume_field-coral_dim$volume_pg #actual overestimation

coral_dim$prop_est <- (coral_dim$volume_field - coral_dim$volume_pg)/coral_dim$volume_pg #proportion overestimated 

coral_dim <- coral_dim %>% 
  dplyr::select(coral_id, branch, cafi, size_class, volume_field, volume_pg, SA, height_pg, length_pg, width_pg, avg_w_cm, interstitial_space, SAV, convexity, packing, sphericity, diff_est, prop_est) %>%
  filter(cafi=="empty")
  
#coral_dim$branch <- ifelse(coral_dim$convexity>=0.5, "tight","wide") #classifies wide and tight branching coral based on convexity measurement

Merge invertebrate data with coral morphological data.

Dataframes were combined and seperate columns created for log and square root transformations.

#Create dataframe called cafi_coral and merge coral and cafi data
cafi_coral <- merge(cafi_summarized2, coral_dim, by ="coral_id") %>%
  drop_na(interstitial_space) #%>%
  #dplyr::select(coral_id, branch, cafi, size_class, volume_field, volume_pg, SA, height_pg, length_pg, width_pg, avg_w_cm, interstitial_space, SAV, convexity, packing, sphericity, diff_est, prop_est, num_cafi, cafi_present, cafi_richness, sw) 

cafi_coral$branch <- ifelse(cafi_coral$convexity>=0.5, "tight","wide") #classifies wide and tight branching coral based on convexity measurement

#Log transform variables
cafi_coral$log_num_cafi <- log(cafi_coral$num_cafi)
cafi_coral$log_volume_field <- log(cafi_coral$volume_field)
cafi_coral$log_volume_pg <- log(cafi_coral$volume_pg)
cafi_coral$log_avg_w <- log(cafi_coral$avg_w_cm)
cafi_coral$log_interstitial_space <- log(cafi_coral$interstitial_space)
cafi_coral$log_SA <- log(cafi_coral$SA)
cafi_coral$log_cafi_richness <- log(cafi_coral$cafi_richness)
cafi_coral$log_sphericity <- log(cafi_coral$sphericity)

#Sqrt transform variables
cafi_coral$sqrt_num_cafi <- sqrt(cafi_coral$num_cafi)
cafi_coral$sqrt_volume_pg <- sqrt(cafi_coral$volume_pg)
cafi_coral$sqrt_avg_w_cm <- sqrt(cafi_coral$avg_w_cm)
cafi_coral$sqrt_volume_field <- sqrt(cafi_coral$volume_field)
cafi_coral$sqrt_interstitial_space <- sqrt(cafi_coral$interstitial_space)
cafi_coral$sqrt_SA <- sqrt(cafi_coral$SA)
cafi_coral$sqrt_cafi_richness <- sqrt(cafi_coral$cafi_richness)
cafi_coral$sqrt_sphericity <- sqrt(cafi_coral$sphericity)

Histogram

Look for unimodal or bimodal curve in data for average branch distance, interstitial space and convexity. Compare these curves to branch distances defined in the field to branch distances defined by convexity measuremnt of 0.5

#Load plyr package
library(plyr)

## ------------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## ------------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following object is masked from 'package:faraway':
## 
##     ozone

## The following object is masked from 'package:ggpubr':
## 
##     mutate

## The following objects are masked from 'package:reshape':
## 
##     rename, round_any

## The following object is masked from 'package:purrr':
## 
##     compact

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:here':
## 
##     here

#Classification based on convexity measuremnts
mu <- ddply(cafi_coral, "branch", 
            summarise, grp.mean=mean(avg_w_cm, na.rm=TRUE))
d1 <- ggplot(cafi_coral, aes(x=avg_w_cm, fill=branch))+
  geom_density(alpha=0.4)+
  geom_vline(data=mu, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  labs(x="Branch distance")+
  theme(legend.position = "top")

mu2 <- ddply(cafi_coral, "branch",
             summarise, grp.mean=mean(interstitial_space, na.rm=TRUE))
d2 <- ggplot(cafi_coral, aes(x=interstitial_space, fill=branch))+
  geom_density(alpha=0.4)+
  geom_vline(data=mu2, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  labs(x="Interstitial space")+
  theme(legend.position = "top")

mu3 <- ddply(cafi_coral, "branch",
             summarise, grp.mean=mean(convexity, na.rm=TRUE))
d3 <- ggplot(cafi_coral, aes(x=convexity, fill=branch))+
  labs(x="Convexity")+
  geom_density(alpha=0.4)+
  geom_vline(data=mu3, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  theme(legend.position = "top")

#Classification based on qualitative analysis in the field
mu4 <- ddply(coral_dim, "branch",
             summarise, grp.mean=mean(avg_w_cm, na.rm=TRUE))
d4 <- ggplot(coral_dim, aes(x=avg_w_cm, fill=branch))+
  geom_density(alpha=0.4)+
  geom_vline(data=mu4, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  labs(x="Branch distance")+
  theme(legend.position = "top")

mu5 <- ddply(coral_dim, "branch",
             summarise, grp.mean=mean(interstitial_space, na.rm=TRUE))
d5 <- ggplot(coral_dim, aes(x=interstitial_space, fill=branch))+
  geom_density(alpha=0.4)+
  geom_vline(data=mu5, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  labs(x="Interstitial space")+
  theme(legend.position = "top")

mu6 <- ddply(coral_dim, "branch",
             summarise, grp.mean=mean(convexity, na.rm=TRUE))
d6 <- ggplot(coral_dim, aes(x=convexity, fill=branch))+
  geom_density(alpha=0.4)+
  geom_vline(data=mu6, aes(xintercept=grp.mean, color=branch),
             linetype="dashed")+
  labs(x="Convexity")+
  theme(legend.position = "top")

fig1 <- ggarrange(d1,d2,d3,d4,d5,d6,
          ncol = 3, nrow = 2)

annotate_figure(fig1,
                fig.lab = "Figure 1", fig.lab.face = "bold") 

Visualize correlations

Non-transformed and transformed morphological measurements correlated using *pairs.panels().

#Visualize correlations between variables
cor_var <- cafi_coral %>% 
  dplyr::select(num_cafi,volume_pg,height_pg,volume_field,SA,avg_w_cm,interstitial_space,sphericity,convexity,packing) 
pairs.panels(cor_var, 
             method="pearson", 
             hist.col = "white",
             density=TRUE,
             ellipses=FALSE,
             scale=TRUE)

#Use sqrt transformed data
cor_var2 <- cafi_coral %>%
  dplyr::select(sqrt_num_cafi,sqrt_volume_pg,sqrt_volume_field,sqrt_SA,sqrt_avg_w_cm,sqrt_interstitial_space,sqrt_sphericity,convexity,packing)
pairs.panels(cor_var2, 
             method="pearson", 
             hist.col = "white",
             density=TRUE,
             ellipses=FALSE,
             scale=TRUE)

#Use log transformed data
cor_var3 <- cafi_coral %>%
  dplyr::select(log_num_cafi,log_volume_pg,log_volume_field,log_SA,log_avg_w,log_interstitial_space,log_sphericity,convexity,packing)
pairs.panels(cor_var3, 
             method="pearson", 
             hist.col = "white",
             density=TRUE,
             ellipses=FALSE,
             scale=TRUE)

Log transformed field and software estimated volume and log transformed software estimated SA and interstitial space are all highly correlated to eachother. Packing is also highly correlated to volume and SA.

Paired dependent t-test

Here, we want to compare the differences between our software and manual measurements of volume, height, length and width. Manual measurements for coral volume were estimated from an elipsoid.

dim_pg <- coral_pg %>% 
  dplyr::select(coral_id, height_pg, length_pg, width_pg, volume_pg) %>%
  mutate(height_pg = height_pg*100,
         length_pg = length_pg*100,
         width_pg = width_pg*100,
         volume_pg = volume_pg*10^6) 

dim_field <- coral_field %>%
  dplyr::select(coral_id, height_field, length_field, width_field, volume_field) 

dim_bind <- merge(dim_pg, dim_field, by="coral_id") %>%
  drop_na() %>%
  summarise(coral_id,volume_pg,volume_field,
            height_pg,height_field,
            length_pg,length_field,
            width_pg,width_field, 
            prop_est=(volume_field-volume_pg)/volume_pg, #proportion overestimation of elipsoide compared to software volume
            diff_est=volume_field-volume_pg) #actual overestimation difference 

Since the corals are not independent, that is, they were taken on the same coral but with different methods, we will consider a paired t-test with dependent samples. Height, length and width measurements are also subjective since they may have been taken from different areas on the coral. We will first check for normallity assumptions. Normality will be checked with a qqplot.

par(mfrow=c(3,3))
qqPlot(dim_bind$height_pg, dist="norm")

## [1]  1 58

qqPlot(dim_bind$height_field, dist="norm")

## [1]  1 56

qqPlot(dim_bind$length_pg, dist="norm")

## [1]  6 57

qqPlot(dim_bind$length_field, dist="norm")

## [1]  6 57

qqPlot(dim_bind$width_pg, dist="norm")

## [1] 19  1

qqPlot(dim_bind$width_field, dist="norm")

## [1]  1 23

par(mfrow=c(2,2))

par(mfrow=c(2,2))
qqPlot(dim_bind$volume_pg, dist="norm")

## [1] 58 23

qqPlot(dim_bind$volume_field, dist="norm")

## [1]  1 36

par(mfrow=c(1,1))

#Volume data is not normal, try sqrt transformation
dim_bind$sqrt_volume_pg <- sqrt(dim_bind$volume_pg)
dim_bind$sqrt_volume_field <- sqrt(dim_bind$volume_field)

par(mfrow=c(2,2))
qqPlot(dim_bind$sqrt_volume_pg, dist="norm")

## [1] 58 23

qqPlot(dim_bind$sqrt_volume_field, dist="norm")

## [1]  1 36

par(mfrow=c(1,1))

#sqrt transformed volume for both measurements is normal, move to paired t-test

t.test(dim_bind$sqrt_volume_pg,dim_bind$sqrt_volume_field, paired=TRUE)

## 
##  Paired t-test
## 
## data:  dim_bind$sqrt_volume_pg and dim_bind$sqrt_volume_field
## t = -8.8976, df = 58, p-value = 1.957e-12
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -29.08742 -18.40334
## sample estimates:
## mean of the differences 
##               -23.74538

#ddply(dim_bind, .(sqrt_volume_pg,sqrt_volume_field), summarize,
t.test(dim_bind$height_pg,dim_bind$height_field, paired=TRUE)

## 
##  Paired t-test
## 
## data:  dim_bind$height_pg and dim_bind$height_field
## t = -5.6455, df = 58, p-value = 5.194e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.483686 -1.183436
## sample estimates:
## mean of the differences 
##               -1.833561

t.test(dim_bind$length_pg,dim_bind$length_field, paired=TRUE)

## 
##  Paired t-test
## 
## data:  dim_bind$length_pg and dim_bind$length_field
## t = -0.0097688, df = 58, p-value = 0.9922
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.5147722  0.5097722
## sample estimates:
## mean of the differences 
##                 -0.0025

t.test(dim_bind$width_pg,dim_bind$width_field, paired=TRUE)

## 
##  Paired t-test
## 
## data:  dim_bind$width_pg and dim_bind$width_field
## t = -4.2154, df = 58, p-value = 8.836e-05
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.4710549 -0.5237823
## sample estimates:
## mean of the differences 
##              -0.9974186

We reject the null hypothesis that the true difference in means is equal to zero for all measurements except for length.

Comparing morphological measurements between tight and wide branching corals

#Log SA
g1 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_SA, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="none")+
  labs(y=expression(paste("Log SA cm"^{2})))

#Log avg. width
g2 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_avg_w, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="none")+
  labs(y="Log branch distance cm")

#Log IS
g3 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_interstitial_space, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="none")+
  labs(y=expression(paste("Log space cm"^{3})))

#Log Sphericity
g4 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_sphericity, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="none")+
  labs(y="Log sphericity")

#Convexity
g5 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=convexity, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="none")+
  labs(y="Convexity", x=expression(paste("Log Vol. cm"^{3})))

#Packing
g6 <- ggplot(cafi_coral, aes(x=log_volume_pg, y=packing, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_classic()+
  theme(axis.title.x=element_blank(),
        legend.position="bottom")+
  labs(y="Packing")

#extract legend
#https://github.com/hadley/ggplot2/wiki/Share-a-legend-between-two-ggplot2-graphs
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(g6)
       
fig2 <- ggarrange(g1,g2,g3,g4,g5,g6,
                     ncol = 3, nrow=2)
  
annotate_figure(fig2,
                top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
                bottom = text_grob(expression(paste("Log Vol. cm"^{3})), color = "black",
                                   face = "bold", size = 12),
                fig.lab = "Figure 2", fig.lab.face = "bold")  

Graphically analyze software estimated volume and elipsoid volume

The below graphs compare the software estimated volume, height, length and width measurements to the same coral’s manual measurements. A one to one line was fit to show how differences in measuremets vary with increasing coral size.

#Volume
#summary(lm(formula = volume_field ~ volume_pg, data = dim_bind))
volume <- ggplot(dim_bind, aes(x=volume_pg, y=volume_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  scale_y_continuous(breaks=seq(0,max(dim_bind$volume_field), by = 5000))+
  scale_x_continuous(breaks=seq(0,max(dim_bind$volume_pg), by = 2000))+
  theme_bw()+
  labs(y= expression(paste("elipsoid V cm"^{3})), x= expression(paste("software V cm"^{3})))+
  geom_abline(intercept = 756.3176, slope = 1, color='black') #x-intercept pulled from summary data

#Height
#summary(lm(formula = height_field ~ height_pg, data = dim_bind))
height <- ggplot(dim_bind, aes(x=height_pg, y=height_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y= "Manual h cm", x= "Software h cm")+
  geom_abline(intercept = 2.44284, slope = 1, color='black') 

#Length
#summary(lm(formula = length_field ~ length_pg, data = dim_bind))
length <- ggplot(dim_bind, aes(x=length_pg, y=length_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y= "Manual l cm", x= "Software l cm")+
  geom_abline(intercept = 0.52596, slope = 1, color='black')

#Width
#summary(lm(formula = width_field ~ width_pg, data = dim_bind))
width <- ggplot(dim_bind, aes(x=width_pg, y=width_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y= "Manual w cm", x= "Software w cm")+
  geom_abline(intercept = 2.13398, slope = 1, color='black')
  
fig3 <- ggarrange(volume, height, length, width,
          ncol = 2, nrow = 2)

annotate_figure(fig3,
                top = text_grob("Manual v. software measurements", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 3", fig.lab.face = "bold")  

The below graphs show the proportion and actual differences of software volume from that of elipsoid volume. This change is shown with increasing software coral volume which we will assume is the closest to the actual volume of the coral.

#Actual overestimation
plot_est <- ggplot(cafi_coral, aes(x=volume_pg, y=diff_est, shape=branch, color=branch, linetype=branch))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_bw()+
   labs(y= expression(paste("elipsoid overestim. cm"^{3})), x= expression(paste("software V cm"^{3})))+
  theme(legend.position="none")

#Linear Model
lm_plot_est <- lm(diff_est~volume_pg+branch, data=cafi_coral) 

#Check normality, homoscedasticity, linearity
par(mfrow=c(2,2))
plot(lm_plot_est)

par(mfrow=c(1,1))

#Proportional overestimation
plot_prop <- ggplot(cafi_coral, aes(x=volume_pg, y=prop_est, shape=branch, color=branch, linetype=branch))+
  geom_point()+
  geom_smooth(method=lm, se=TRUE, aes(fill=branch))+
  theme_bw()+
   labs(y= "Proportion overestim.", x= expression(paste("software V cm"^{3})))

#Linear Model
lm_plot_prop <- lm(prop_est~volume_pg+branch, data=cafi_coral) 

#Check normality, homoscedasticity, linearity
par(mfrow=c(2,2))
plot(lm_plot_prop)

par(mfrow=c(1,1))

fig1 <- ggarrange(plot_est, plot_prop,
          ncol = 2, nrow = 1)

annotate_figure(fig1,
                top = text_grob("Actual vs. proportional estimates", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 1", fig.lab.face = "bold")  

ggsave("fig1.pdf", plot = fig1, path="figs", width=7, height=4)

#Create table with the linear model parameters
pl <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software V")

pl2 <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software V")

tab_model(lm_plot_est,lm_plot_prop,
          pred.labels = c(pl,pl2),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          dv.labels = c("Overestim.", "Proportion over-est."))

	Overestim.				Proportion over-est.
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	-2587.16	2075.50	-6853.42 – 1679.09	0.224	0.51	0.27	-0.05 – 1.07	0.075
Software V	1.14	0.26	0.61 – 1.67	<0.001	0.00	0.00	-0.00 – 0.00	0.333
branchwide	4843.95	2066.56	596.07 – 9091.82	0.027	1.06	0.27	0.50 – 1.62	0.001
Observations	29				29
R2 / R2 adjusted	0.434 / 0.390				0.371 / 0.323

Measuring the volume of an elipsoid overestimates actual coral volume approximately by 50% across all size ranges and overestimates this volume by a difference of approximately \(4280.195 cm^3\) Looking at the plots, as our best estimate of coral volume increases, the difference between the elipsoid volume and our best estimate also increases in a positive relationship. From our second plot of overestimation proportion with increasing coral volume as the predictor variable, there is no evidence of the proportion or \(software V/elipsoid V\) increasing with increasing best estimate volume with a \(p-value=0.061\). The proportion stayed constant at approximately 50% overestimation of our best estimate of the actual coral volume. ########################################################################################################################

The below graphs compare the sqrt and log transformed software estimated volume to the same coral’s sqrt and log transformed elipsoid volume. A one to one line was fit to show how differences in measuremets vary with increasing coral size.

#Sqrt volume
#summary(lm(formula = sqrt_volume_field ~ sqrt_volume_pg, data = cafi_coral))
sqrt_volume <- ggplot(cafi_coral, aes(x=sqrt_volume_pg, y=sqrt_volume_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y= expression(paste("sqrt elipsoid V cm"^{3})), x= expression(paste("sqrt software V cm"^{3})))+
  geom_abline(intercept =  4.207, slope = 1, color='blue') #x-intercept pulled from summary data

#summary(lm(formula = log_volume_field ~ log_volume_pg, data = cafi_coral))
log_volume <- ggplot(cafi_coral, aes(x=log_volume_pg, y=log_volume_field))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y= expression(paste("log elipsoid V cm"^{3})), x= expression(paste("log software V cm"^{3})))+
  geom_abline(intercept =   0.80673, slope = 1, color='blue') #x-intercept pulled from summary data

fig4 <- ggarrange(volume, sqrt_volume, log_volume,
          ncol = 3, nrow = 1)

annotate_figure(fig4,
                top = text_grob("Compare transformations", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 4", fig.lab.face = "bold")  

The below graphs represent how invertebrate abundance and richness relates to both volume measurements. We will test how well our software estimated method of calculating volume and interstitial space relates to invertebrate counts. We will first assess \(p<0.05\) to test if the linear regression slope was significantly different from zero. Correlation of dependent and independent variables \(R^2\) was also included in the tables.

#Plot invert abundance to elipsoid volume
invert_field <- ggplot(cafi_coral, aes(x=volume_field, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  labs(y="invert abundance", x= expression(paste("elipsoid V cm"^{3})))

#Check qqplots for normality 
lm_invert_field <- lm(num_cafi~volume_field, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_invert_field)

par(mfrow=c(1,1))

#Plot invert abundance to software estimated volume
invert_pg <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  labs(y="invert abundance", x= expression(paste("software V cm"^{3})))

#Check qqplots for normality 
lm_invert_pg <- lm(num_cafi~volume_pg, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_invert_pg)

par(mfrow=c(1,1))

#Plot invert richness to software estimated volume
richness_field <- ggplot(cafi_coral, aes(x=volume_field, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  labs(y="invert richness", x= expression(paste("elipsoid V cm"^{3})))

#Check qqplots for normality 
lm_richness_field <- lm(cafi_richness~volume_field, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_richness_field)

par(mfrow=c(1,1))

#Plot invert richness to software estimated volume
richness_pg <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  labs(y="invert richness", x= expression(paste("software V cm"^{3})))

#Check qqplots for normality
lm_richness_pg <- lm(cafi_richness~volume_pg, data=cafi_coral)
par(mfrow=c(2,2))
plot(lm_richness_pg)

par(mfrow=c(1,1))

fig2 <- ggarrange(invert_field, invert_pg, richness_field, richness_pg,
          ncol = 2, nrow = 2)

annotate_figure(fig2,
                top = text_grob("Invert Abundance and Richness", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 2", fig.lab.face = "bold") 

ggsave("fig2.pdf", plot = fig2, path="figs", width=5, height=4)

#HTML table comparing volume methods and abundance
pl3 <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid volume")
pl4 <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software est. volume")

tab_model(lm_invert_field,lm_invert_pg,
          pred.labels = c(pl3,pl4),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Field volume vs. abundance","Software volume vs. abundance"))

	Field volume vs. abundance				Software volume vs. abundance
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	13.69118	3.08747	7.35622 – 20.02615	<0.001	11.54653	2.92781	5.53916 – 17.55390	0.001
Elipsoid volume	0.00052	0.00024	0.00002 – 0.00102	0.043
Software est. volume					0.00169	0.00053	0.00060 – 0.00277	0.004
Observations	29				29
R2 / R2 adjusted	0.143 / 0.111				0.275 / 0.248

#HTML table comparing volume methods and richness
pl5 <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid volume")
pl6 <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software est. volume")

tab_model(lm_richness_field,lm_richness_pg,
          pred.labels = c(pl5,pl6),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Field volume vs. richness","Software volume vs. richness"))

	Field volume vs. richness				Software volume vs. richness
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	6.24931	0.85828	4.48826 – 8.01037	<0.001	5.87648	0.84250	4.14782 – 7.60514	<0.001
Elipsoid volume	0.00019	0.00007	0.00005 – 0.00033	0.011
Software est. volume					0.00051	0.00015	0.00019 – 0.00082	0.003
Observations	29				29
R2 / R2 adjusted	0.218 / 0.189				0.291 / 0.265

When comparing elipsoid volume to software estimated volume for invertebrate richness, we see significant p-values for both as represented in the above tables. Software estimated volume measurements and invert. richness show a lower p-value indicating a greater slope from zero as well as having a higher correlation \(R^2\) value.

Graphically analyze software estimated space and avg. branch distance

Software estimated available space calculated from convex hull and coral volume was related to invertebrate abundance. The same was repeated for the traditional method of average branch distance.

Dependent and independent variables were also compared after being log transformed.

branch_pg <- ggplot(cafi_coral, aes(x=interstitial_space, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="invert abundance", x= expression(paste("IS cm"^{3})))

lm_branch_pg <- lm(num_cafi~interstitial_space, data=cafi_coral)

log_branch_pg <- ggplot(cafi_coral, aes(x=log_interstitial_space, y=log_num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="Log invert abundance", x= expression(paste("Log IS cm"^{3})))

lm_log_branch_pg <- lm(log_num_cafi~log_interstitial_space, data=cafi_coral)

branch_field <- ggplot(cafi_coral, aes(x=avg_w_cm, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="invert abundance", x= "avg_w_cm")

lm_branch_field <- lm(num_cafi~avg_w_cm, data=cafi_coral)

log_branch_field <- ggplot(cafi_coral, aes(x=log_avg_w, y=log_num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="Log invert abundance", x= "Log avg. width")

lm_log_branch_field <- lm(log_num_cafi~log_avg_w, data=cafi_coral)

fig7 <- ggarrange(branch_field, branch_pg, log_branch_field, log_branch_pg,
          ncol = 2, nrow = 2)

annotate_figure(fig7,
                top = text_grob("Invert abundance and habitat space", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 7", fig.lab.face = "bold")  

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  interstitial_space = "IS")

tab_model(lm_branch_pg,
          pred.labels = pl,
          dv.labels = "Invertebrate abundance")

	Invertebrate abundance
Predictors	Estimates	CI	p
Intercept	16.01	9.85 – 22.17	<0.001
IS	0.00	-0.00 – 0.00	0.250
Observations	29
R2 / R2 adjusted	0.049 / 0.013

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  log_interstitial_space = "Log IS")

tab_model(lm_log_branch_pg,
          pred.labels = pl,
          dv.labels = "Invertebrate abundance")

	Invertebrate abundance
Predictors	Estimates	CI	p
Intercept	0.56	-0.59 – 1.70	0.329
Log IS	0.28	0.13 – 0.43	0.001
Observations	29
R2 / R2 adjusted	0.354 / 0.330

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  avg_w_cm = "Avg. branch width")

tab_model(lm_branch_field,
          pred.labels = pl,
          dv.labels = "Invertebrate abundance")

	Invertebrate abundance
Predictors	Estimates	CI	p
Intercept	17.06	5.45 – 28.66	0.006
Avg. branch width	0.60	-5.02 – 6.21	0.829
Observations	29
R2 / R2 adjusted	0.002 / -0.035

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  log_avg_w = "Log avg. width")

tab_model(lm_log_branch_field,
          pred.labels = pl,
          dv.labels = "Invertebrate abundance")

	Invertebrate abundance
Predictors	Estimates	CI	p
Intercept	2.54	2.10 – 2.98	<0.001
Log avg. width	0.23	-0.42 – 0.87	0.476
Observations	29
R2 / R2 adjusted	0.019 / -0.017

Log transformed data appeared to not violate normality assumptions so was used in this analysis. Relating interstitial space to invertebrate abundance gave a signficant p-value while average width did not. The difference in p-values and \(R^2\) correlation values were significant. This may provide evidence to average branch width not adequately capturing available space for invertebrates when compared to a 3D software estimated measurement.

branch_rich_pg <- ggplot(cafi_coral, aes(x=interstitial_space, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="invert richness", x= expression(paste("IS cm"^{3})))

lm_branch_rich_pg <- lm(cafi_richness~interstitial_space, data=cafi_coral)

log_branch_rich_pg <- ggplot(cafi_coral, aes(x=log_interstitial_space, y=log_cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="Log invert richness", x= expression(paste("Log IS cm"^{3})))

lm_log_branch_rich_pg <- lm(log_cafi_richness~log_interstitial_space, data=cafi_coral)

branch_rich_field <- ggplot(cafi_coral, aes(x=avg_w_cm, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="invert richness", x= "avg_w_cm")

lm_branch_rich_field <- lm(cafi_richness~avg_w_cm, data=cafi_coral)

log_branch_rich_field <- ggplot(cafi_coral, aes(x=log_avg_w, y=log_cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_bw()+
  labs(y="Log invert richness", x= "Log avg. width")

lm_log_branch_rich_field <- lm(log_cafi_richness~log_avg_w, data=cafi_coral)

fig8 <- ggarrange(branch_rich_field, branch_rich_pg, log_branch_rich_field, log_branch_rich_pg,
          ncol = 2, nrow = 2)

annotate_figure(fig8,
                top = text_grob("Invert richness and habitat space", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 8", fig.lab.face = "bold")  

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  interstitial_space = "IS")

tab_model(lm_branch_rich_pg,
          pred.labels = pl,
          dv.labels = "Invertebrate richness")

	Invertebrate richness
Predictors	Estimates	CI	p
Intercept	6.83	5.11 – 8.54	<0.001
IS	0.00	-0.00 – 0.00	0.053
Observations	29
R2 / R2 adjusted	0.131 / 0.099

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  log_interstitial_space = "Log IS")

tab_model(lm_log_branch_rich_pg,
          pred.labels = pl,
          dv.labels = "Invertebrate richness")

	Invertebrate richness
Predictors	Estimates	CI	p
Intercept	0.33	-0.47 – 1.14	0.403
Log IS	0.21	0.11 – 0.32	<0.001
Observations	29
R2 / R2 adjusted	0.390 / 0.368

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  avg_w_cm = "Avg. branch width")

tab_model(lm_branch_rich_field,
          pred.labels = pl,
          dv.labels = "Invertebrate richness")

	Invertebrate richness
Predictors	Estimates	CI	p
Intercept	6.50	3.17 – 9.83	<0.001
Avg. branch width	0.73	-0.88 – 2.34	0.361
Observations	29
R2 / R2 adjusted	0.031 / -0.005

#HTML table from regression results
pl <- c(
  '(Intercept)' = "Intercept",
  log_avg_w = "Log avg. width")

tab_model(lm_log_branch_rich_field,
          pred.labels = pl,
          dv.labels = "Invertebrate richness")

	Invertebrate richness
Predictors	Estimates	CI	p
Intercept	1.81	1.50 – 2.13	<0.001
Log avg. width	0.25	-0.22 – 0.71	0.285
Observations	29
R2 / R2 adjusted	0.042 / 0.007

Log transformed data appeared to not violate normality assumptions so was used in this analysis. Relating interstitial space to invertebrate richness gave a signficant p-value while average width did not. The difference in p-values and \(R^2\) correlation values were significant. This may provide evidence to average branch width not adequately capturing available space for invertebrates when compared to a 3D software estimated measurement.

Compare mean invertebrate abundance and species richness of tight and wide corals

Visualize data using boxplots.

# Boxplot of invert abundance
log_invert_box <- ggplot(cafi_coral, aes(x= branch, y=log_num_cafi, color = branch))+
  geom_boxplot(width = 0.5)+
  geom_jitter()+ #adds data on top of box plot
  scale_color_manual(values = c("firebrick1", "steelblue2"))+
  scale_y_continuous(limits= c(0, max(cafi_coral$log_num_cafi + 1)))+
  theme_classic()+
  labs(y="Log invert. abundance", x="Coral branching", title = "Log invert. abundance boxplot")

log_invert_box

# Boxplot of invertebrate species richness per coral for both tight and wide branching
richness_box <- ggplot(cafi_coral, aes(x= branch, y=log_cafi_richness, color = branch))+
  geom_boxplot(width = 0.5)+
  geom_jitter()+
  scale_color_manual(values = c("firebrick1", "steelblue2"))+
  scale_y_continuous(limits= c(0, max(cafi_coral$log_cafi_richness + 1)))+
  theme_classic()+
  labs(y="Invert. sp. richness", x="Coral branching", title = "Invert. sp. richness boxplot")

richness_box

# Boxplot of invertebrate diversity for both tight and wide branching coral
diversity_box <- ggplot(cafi_coral, aes(x= branch, y=sw, color = branch))+
  geom_boxplot(width = 0.5)+
  geom_jitter()+
  scale_color_manual(values = c("firebrick1", "steelblue2"))+
  scale_y_continuous(limits= c(0, max(cafi_coral$sw + 0.5)), breaks = c(0, 0.5, 1, 1.5, 2, 2.5, 3))+
  theme_classic()+
  labs(y="Invert. diversity", x="Coral branching", title = "Invert. diversity boxplot")

diversity_box

Data was log transformed except for diversity so as not to violate normality assumptions. It appears that there is no difference in abundance, species richness or diversity.

Before conducting a two sample t-test we checked for normality using qq-plots and Shapiro Wilk test and a Levene’s test for equal variance.

#Q-Q plot separated by group 
with(cafi_coral, qqPlot(num_cafi[branch == "wide"], dist="norm")) 

## [1] 13 14

with(cafi_coral, qqPlot(num_cafi[branch == "tight"], dist="norm")) 

## [1] 12  6

#Shapiro Wilk seperated by group
with(cafi_coral, shapiro.test(num_cafi[branch == "wide"])) 

## 
##  Shapiro-Wilk normality test
## 
## data:  num_cafi[branch == "wide"]
## W = 0.92775, p-value = 0.1995

with(cafi_coral, shapiro.test(num_cafi[branch == "tight"])) 

## 
##  Shapiro-Wilk normality test
## 
## data:  num_cafi[branch == "tight"]
## W = 0.76166, p-value = 0.003547

With a \(p>0.05\) we fail to reject the null hypothesis that our data is normally distributed for wide branching corals. With a \(p<0.05\) we reject the null hypothesis that our data is normally distributed for tight branching corals. Because we are working with smaller data sets, we cannot use the central limit theorem.

#Levene's test: response variable, group variable
leveneTest(cafi_coral$num_cafi, cafi_coral$branch)

## Warning in leveneTest.default(cafi_coral$num_cafi, cafi_coral$branch):
## cafi_coral$branch coerced to factor.

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.2017 0.6569
##       27

With a \(p>0.05\) we fail to reject the null hypothesis that our data is normally distributed for tight and wide branching corals. Both groups of corals that have tight and wide branching morphologies have equal variance.

Because tight branching coral is not normally distributed, we could either perform a transformation or non-parametric test. We will try log transforming our data first. A column log transforming the number of inverts has already been created

#Q-Q plot separated by group 
with(cafi_coral, qqPlot(log_num_cafi[branch == "wide"], dist="norm")) 

## [1] 8 6

with(cafi_coral, qqPlot(log_num_cafi[branch == "tight"], dist="norm")) 

## [1] 12  6

#Shapiro Wilk seperated by group
with(cafi_coral, shapiro.test(log_num_cafi[branch == "wide"])) 

## 
##  Shapiro-Wilk normality test
## 
## data:  log_num_cafi[branch == "wide"]
## W = 0.93506, p-value = 0.2642

with(cafi_coral, shapiro.test(log_num_cafi[branch == "tight"])) 

## 
##  Shapiro-Wilk normality test
## 
## data:  log_num_cafi[branch == "tight"]
## W = 0.88255, p-value = 0.09449

#Levene's test: response variable, group variable
leveneTest(cafi_coral$log_num_cafi, cafi_coral$branch)

## Warning in leveneTest.default(cafi_coral$log_num_cafi, cafi_coral$branch):
## cafi_coral$branch coerced to factor.

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  1  0.4491 0.5084
##       27

# subset data
tight <- subset(cafi_coral, branch == "tight")
wide <- subset(cafi_coral, branch == "wide")

t.test(tight$log_num_cafi, wide$log_num_cafi, var.equal = TRUE)

## 
##  Two Sample t-test
## 
## data:  tight$log_num_cafi and wide$log_num_cafi
## t = 0.9573, df = 27, p-value = 0.3469
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.3005325  0.8262331
## sample estimates:
## mean of x mean of y 
##   2.81695   2.55410

With a \(p>0.05\) we fail to reject the null hypothesis that the difference between the means of the log transformed invert counts are different from zero.

####Compare other photogrametry morphological measurements to invert abundance

#Volume
plot1 <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x=expression(paste("Volume cm"^{3})))

#Sphericity
plot2 <- ggplot(cafi_coral, aes(x=sphericity, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Sphericity")

lm_sphericity <- lm(num_cafi~sphericity, data = cafi_coral)

#Packing
plot3 <- ggplot(cafi_coral, aes(x=packing, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Packing")

lm_packing <- lm(num_cafi~packing, data = cafi_coral)

#Height
plot4 <- ggplot(cafi_coral, aes(x=height_pg, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Height")

lm_height <- lm(num_cafi~height_pg, data = cafi_coral)

#Length
plot5 <- ggplot(cafi_coral, aes(x=length_pg, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Length")

lm_length <- lm(num_cafi~length_pg, data = cafi_coral)

#Width
plot6 <- ggplot(cafi_coral, aes(x=width_pg, y=num_cafi))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Width")

lm_width <- lm(num_cafi~width_pg, data = cafi_coral)

#Figure 3
fig3 <- ggarrange(plot1,plot2,plot3,plot4,plot5,plot6,
                     ncol = 3, nrow=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

annotate_figure(fig3,
                top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
                left = text_grob("Invert. Abundance", color = "black",
                                   rot = 90, face = "bold", size = 12),
                fig.lab = "Figure 3", fig.lab.face = "bold") 

ggsave("fig3.pdf", plot = fig3, path="figs", width=6, height=4)

#Data tables for volume, sphericity, packing
pl7 <- c('(Intercept)' = "Intercept",
  volume_pg = "Software Volume")

pl8 <- c('(Intercept)' = "Intercept",
  sphericity = "Sphericity")

pl9 <- c('(Intercept)' = "Intercept",
  packing = "Packing")

tab_model(lm_invert_pg,lm_sphericity, lm_packing,
          pred.labels = c(pl7,pl8,pl9),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Abundance vs. Software Volume","Abundance vs. Sphericity","Abundance vs. Packing"))

	Abundance vs. Software Volume				Abundance vs. Sphericity				Abundance vs. Packing
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	11.54653	2.92781	5.53916 – 17.55390	0.001	33.93199	9.99889	13.41597 – 54.44801	0.002	-18.68154	12.72981	-44.80096 – 7.43788	0.154
Software Volume	0.00169	0.00053	0.00060 – 0.00277	0.004
Sphericity					-51.45528	31.75640	-116.61402 – 13.70347	0.117
Packing									19.19597	6.53838	5.78032 – 32.61162	0.007
Observations	29				29				29
R2 / R2 adjusted	0.275 / 0.248				0.089 / 0.055				0.242 / 0.214

#Data tables for height, length, width
pl10 <- c('(Intercept)' = "Intercept",
  height_pg = "Height")

pl11 <- c('(Intercept)' = "Intercept",
  length_pg = "Length")

pl12 <- c('(Intercept)' = "Intercept",
  width_pg = "Width")

tab_model(lm_height,lm_length, lm_width,
          pred.labels = c(pl10,pl11,pl12),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Abundance vs. Height","Abundance vs. Length","Abundance vs. Width"))

	Abundance vs. Height				Abundance vs. Length				Abundance vs. Width
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	2.04637	4.54287	-7.27482 – 11.36757	0.656	3.52788	5.19889	-7.13936 – 14.19512	0.503	3.93804	5.02293	-6.36815 – 14.24423	0.440
Height	0.97632	0.24875	0.46593 – 1.48671	0.001
Length					0.56273	0.18292	0.18742 – 0.93804	0.005
Width									0.66226	0.21264	0.22596 – 1.09856	0.004
Observations	29				29				29
R2 / R2 adjusted	0.363 / 0.340				0.260 / 0.232				0.264 / 0.237

ggsave("fig3.pdf", plot = fig3, path="figs", width=6, height=4)

####Compare other photogrametry morphological measurements to invert richness

#Volume
plot7 <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x=expression(paste("Volume cm"^{3})))

#Sphericity
plot8 <- ggplot(cafi_coral, aes(x=sphericity, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Sphericity")

lm_sph_rich <- lm(cafi_richness~sphericity, data=cafi_coral)

#Packing
plot9 <- ggplot(cafi_coral, aes(x=packing, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Packing")

lm_pack_rich <- lm(cafi_richness~packing, data=cafi_coral)

#Height
plot10 <- ggplot(cafi_coral, aes(x=height_pg, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Height")

lm_height_rich <- lm(cafi_richness~height_pg, data=cafi_coral)

#Length
plot11 <- ggplot(cafi_coral, aes(x=length_pg, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Length")

lm_length_rich <- lm(cafi_richness~length_pg, data=cafi_coral)

#Width
plot12 <- ggplot(cafi_coral, aes(x=width_pg, y=cafi_richness))+ 
  geom_point()+
  geom_smooth(method=lm, se=TRUE)+
  theme_classic()+
  theme(axis.title.y=element_blank(),
        legend.position="none")+
  labs(x="Width")

lm_width_rich <- lm(cafi_richness~width_pg, data=cafi_coral)

#Figure 4
fig4 <- ggarrange(plot7,plot8,plot9,plot10,plot11,plot12,
                     ncol = 3, nrow=2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

annotate_figure(fig4,
                top = text_grob("Morphological Measurements", color = "black", face = "bold", size = 14),
                left = text_grob("Invert. Richness", color = "black",
                                   rot = 90, face = "bold", size = 12),
                fig.lab = "Figure 4", fig.lab.face = "bold")

ggsave("fig4.pdf", plot = fig4, path="figs", width=6, height=4)

#Data tables for volume, sphericity, packing
pl7 <- c('(Intercept)' = "Intercept",
  volume_pg = "Software Volume")

pl8 <- c('(Intercept)' = "Intercept",
  sphericity = "Sphericity")

pl9 <- c('(Intercept)' = "Intercept",
  packing = "Packing")

tab_model(lm_richness_pg,lm_sphericity, lm_packing,
          pred.labels = c(pl7,pl8,pl9),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Richness vs. Software Volume","Richness vs. Sphericity","Richness vs. Packing"))

	Richness vs. Software Volume				Richness vs. Sphericity				Richness vs. Packing
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	5.87648	0.84250	4.14782 – 7.60514	<0.001	33.93199	9.99889	13.41597 – 54.44801	0.002	-18.68154	12.72981	-44.80096 – 7.43788	0.154
Software Volume	0.00051	0.00015	0.00019 – 0.00082	0.003
Sphericity					-51.45528	31.75640	-116.61402 – 13.70347	0.117
Packing									19.19597	6.53838	5.78032 – 32.61162	0.007
Observations	29				29				29
R2 / R2 adjusted	0.291 / 0.265				0.089 / 0.055				0.242 / 0.214

#Data tables for height, length, width
pl10 <- c('(Intercept)' = "Intercept",
  height_pg = "Height")

pl11 <- c('(Intercept)' = "Intercept",
  length_pg = "Length")

pl12 <- c('(Intercept)' = "Intercept",
  width_pg = "Width")

tab_model(lm_height_rich,lm_length_rich, lm_width_rich,
          pred.labels = c(pl10,pl11,pl12),
          show.se = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Richness vs. Height","Richness vs. Length","Richness vs. Width"))

	Richness vs. Height				Richness vs. Length				Richness vs. Width
Predictors	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value	Estimates	std. Error	Conf. Int (95%)	P-Value
Intercept	2.70318	1.24142	0.15600 – 5.25037	0.038	3.05376	1.43923	0.10070 – 6.00681	0.043	3.25370	1.39842	0.38439 – 6.12301	0.028
Height	0.31233	0.06797	0.17286 – 0.45181	<0.001
Length					0.18476	0.05064	0.08087 – 0.28866	0.001
Width									0.21441	0.05920	0.09294 – 0.33587	0.001
Observations	29				29				29
R2 / R2 adjusted	0.439 / 0.418				0.330 / 0.305				0.327 / 0.302

Poisson and Negative Binomial GLM

We will now fit a poisson general linear model to our invertebrate count data

poiss_field <- glm(num_cafi~volume_field, family=poisson, data=cafi_coral) #Poisson regression for invert count data using elipsoide volume
summary(poiss_field)

## 
## Call:
## glm(formula = num_cafi ~ volume_field, family = poisson, data = cafi_coral)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -4.2927  -1.8563  -0.4155   1.3963   5.3379  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.671e+00  6.145e-02  43.465  < 2e-16 ***
## volume_field 2.340e-05  3.912e-06   5.982 2.21e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 229.18  on 28  degrees of freedom
## Residual deviance: 197.50  on 27  degrees of freedom
## AIC: 332.46
## 
## Number of Fisher Scoring iterations: 5

poiss_pg <- glm(num_cafi~volume_pg, family=poisson, data=cafi_coral) #Poisson regression for invert count data using software est. volume
summary(poiss_pg)

## 
## Call:
## glm(formula = num_cafi ~ volume_pg, family = poisson, data = cafi_coral)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.5429  -2.1822  -0.2023   1.2582   5.1645  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 2.544e+00  6.653e-02  38.233   <2e-16 ***
## volume_pg   7.738e-05  9.381e-06   8.248   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 229.18  on 28  degrees of freedom
## Residual deviance: 167.64  on 27  degrees of freedom
## AIC: 302.61
## 
## Number of Fisher Scoring iterations: 5

The residual deviance shows overdispersion for both models. Consequently, we will run a negative binomial model to account for overdispersion. To test if elispoid and software estimated models are different, we will perform a Vuong test for non-nested models.

glm.nb1 <- glm.nb(num_cafi~volume_field, data=cafi_coral) #negative binomial model for elipsoid volume and abundance to account for overdispersion
summary(glm.nb1)

## 
## Call:
## glm.nb(formula = num_cafi ~ volume_field, data = cafi_coral, 
##     init.theta = 3.248837444, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.88272  -0.76793   0.00909   0.44732   1.80535  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.548e+00  1.555e-01  16.383  < 2e-16 ***
## volume_field 3.594e-05  1.185e-05   3.033  0.00242 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.2488) family taken to be 1)
## 
##     Null deviance: 36.215  on 28  degrees of freedom
## Residual deviance: 29.666  on 27  degrees of freedom
## AIC: 217.94
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.249 
##           Std. Err.:  0.980 
## 
##  2 x log-likelihood:  -211.938

glm.nb2 <- glm.nb(num_cafi~volume_pg, data=cafi_coral) #negative binomial model for software estimated volume and abundance to account for overdispersion
summary(glm.nb2)

## 
## Call:
## glm.nb(formula = num_cafi ~ volume_pg, data = cafi_coral, init.theta = 3.741397829, 
##     link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.18550  -0.97664  -0.05352   0.47083   1.88660  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 2.489e+00  1.525e-01  16.325  < 2e-16 ***
## volume_pg   8.935e-05  2.635e-05   3.391 0.000696 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(3.7414) family taken to be 1)
## 
##     Null deviance: 40.591  on 28  degrees of freedom
## Residual deviance: 29.749  on 27  degrees of freedom
## AIC: 214.64
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  3.74 
##           Std. Err.:  1.18 
## 
##  2 x log-likelihood:  -208.637

#Data table of results
pl_nb <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid volume")
pl2_nb <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software est. volume")

tab_model(glm.nb1,glm.nb2,
          pred.labels = c(pl_nb,pl2_nb),
          show.se = TRUE,
          show.aicc = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Negative binomial abundance"))

	Negative binomial abundance
Predictors	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value
Intercept	12.78321	0.15554	9.12736 – 18.08719	<0.001	12.05323	0.15249	8.85505 – 16.55706	<0.001
Elipsoid volume	1.00004	0.00001	1.00001 – 1.00007	0.002
Software est. volume					1.00009	0.00003	1.00003 – 1.00015	0.001
Observations	29				29
R2 Nagelkerke	0.283				0.414
AIC	218.898				215.597

#Vuongtest for abundance negative binomial models
library(nonnest2)

## This is nonnest2 0.5-5.
## nonnest2 has not been tested with all combinations of model classes.

vuongtest(glm.nb1,glm.nb2) #compare negative binomial models

## 
## Model 1 
##  Class: negbin 
##  Call: glm.nb(formula = num_cafi ~ volume_field, data = cafi_coral, ...
## 
## Model 2 
##  Class: negbin 
##  Call: glm.nb(formula = num_cafi ~ volume_pg, data = cafi_coral, init.theta = 3.741397829, ...
## 
## Variance test 
##   H0: Model 1 and Model 2 are indistinguishable 
##   H1: Model 1 and Model 2 are distinguishable 
##     w2 = 0.138,   p = 0.108
## 
## Non-nested likelihood ratio test 
##   H0: Model fits are equal for the focal population 
##   H1A: Model 1 fits better than Model 2 
##     z = -0.827,   p = 0.796
##   H1B: Model 2 fits better than Model 1 
##     z = -0.827,   p = 0.2042

We will now fit a poisson general linear model to our invertebrate richness data

poiss_field_r <- glm(cafi_richness~volume_field, family=poisson, data=cafi_coral) #Poisson regression for richness using elipsoide volume
summary(poiss_field_r)

## 
## Call:
## glm(formula = cafi_richness ~ volume_field, family = poisson, 
##     data = cafi_coral)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.14226  -0.96518  -0.04043   0.72953   2.43666  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.869e+00  9.294e-02  20.116  < 2e-16 ***
## volume_field 2.004e-05  6.110e-06   3.279  0.00104 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 49.269  on 28  degrees of freedom
## Residual deviance: 39.591  on 27  degrees of freedom
## AIC: 154
## 
## Number of Fisher Scoring iterations: 4

poiss_pg_r <- glm(cafi_richness~volume_pg, family=poisson, data=cafi_coral) #Poisson regression for invert richness using software est. volume
summary(poiss_pg_r)

## 
## Call:
## glm(formula = cafi_richness ~ volume_pg, family = poisson, data = cafi_coral)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.2600  -0.9414   0.1766   0.5776   2.4931  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.815e+00  9.864e-02  18.396  < 2e-16 ***
## volume_pg   5.619e-05  1.486e-05   3.783 0.000155 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 49.269  on 28  degrees of freedom
## Residual deviance: 36.051  on 27  degrees of freedom
## AIC: 150.46
## 
## Number of Fisher Scoring iterations: 4

The residual deviance shows overdispersion for both models. Consequently, we will run a negative binomial model to account for overdispersion. To test if elispoid and software estimated models are different, we will perform a Vuong test for non-nested models.

glm.nb3 <- glm.nb(cafi_richness~volume_field, data=cafi_coral) #negative binomial model for elipsoid volume and richness to account for overdispersion
summary(glm.nb3)

## 
## Call:
## glm.nb(formula = cafi_richness ~ volume_field, data = cafi_coral, 
##     init.theta = 20.91987063, link = log)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.91079  -0.83442  -0.01783   0.60680   2.03014  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.855e+00  1.090e-01  17.025  < 2e-16 ***
## volume_field 2.152e-05  7.547e-06   2.852  0.00435 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(20.9199) family taken to be 1)
## 
##     Null deviance: 35.986  on 28  degrees of freedom
## Residual deviance: 28.700  on 27  degrees of freedom
## AIC: 154.19
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  20.9 
##           Std. Err.:  19.5 
## 
##  2 x log-likelihood:  -148.194

glm.nb4 <- glm.nb(cafi_richness~volume_pg, data=cafi_coral) #negative binomial model for software estimated volume and richness to account for overdispersion
summary(glm.nb4)

## 
## Call:
## glm.nb(formula = cafi_richness ~ volume_pg, data = cafi_coral, 
##     init.theta = 33.32422167, link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9972  -0.8639   0.1459   0.5140   2.2036  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 1.810e+00  1.090e-01  16.602  < 2e-16 ***
## volume_pg   5.725e-05  1.703e-05   3.361 0.000776 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(33.3242) family taken to be 1)
## 
##     Null deviance: 39.936  on 28  degrees of freedom
## Residual deviance: 29.257  on 27  degrees of freedom
## AIC: 151.72
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  33.3 
##           Std. Err.:  45.3 
## 
##  2 x log-likelihood:  -145.723

#Data table of results
pl3_nb <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid volume")
pl4_nb <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software est. volume")

tab_model(glm.nb3,glm.nb4,
          pred.labels = c(pl3_nb,pl4_nb),
          show.se = TRUE,
          show.aicc = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Negative binomial richness"))

	Negative binomial richness
Predictors	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value
Intercept	6.39297	0.10897	5.13025 – 7.92216	<0.001	6.10986	0.10902	4.91180 – 7.54485	<0.001
Elipsoid volume	1.00002	0.00001	1.00001 – 1.00004	0.004
Software est. volume					1.00006	0.00002	1.00002 – 1.00009	0.001
Observations	29				29
R2 Nagelkerke	0.313				0.412
AIC	155.154				152.683

#Vuong test for richness negative binomial models
vuongtest(glm.nb3,glm.nb4)

## 
## Model 1 
##  Class: negbin 
##  Call: glm.nb(formula = cafi_richness ~ volume_field, data = cafi_coral, ...
## 
## Model 2 
##  Class: negbin 
##  Call: glm.nb(formula = cafi_richness ~ volume_pg, data = cafi_coral, ...
## 
## Variance test 
##   H0: Model 1 and Model 2 are indistinguishable 
##   H1: Model 1 and Model 2 are distinguishable 
##     w2 = 0.117,   p = 0.233
## 
## Non-nested likelihood ratio test 
##   H0: Model fits are equal for the focal population 
##   H1A: Model 1 fits better than Model 2 
##     z = -0.670,   p = 0.748
##   H1B: Model 2 fits better than Model 1 
##     z = -0.670,   p = 0.2515

#Plot GLM for abundance
plot_nb1 <- ggplot(cafi_coral, aes(x=volume_field, y=num_cafi))+
  geom_point()+
  theme_classic()+
  geom_smooth(method="glm.nb", se=TRUE)+
  labs(y="Invert. Abundance", x= expression(paste("Elipsoid Volume cm"^{3})))


plot_nb2 <- ggplot(cafi_coral, aes(x=volume_pg, y=num_cafi))+
  geom_point()+
  theme_classic()+
  geom_smooth(method="glm.nb", se=TRUE)+
  labs(y="Invert. Abundance", x= expression(paste("Software Volume cm"^{3})))

#Plot GLM for richness
plot_nb3 <- ggplot(cafi_coral, aes(x=volume_field, y=cafi_richness))+
  geom_point()+
  theme_classic()+
  geom_smooth(method="glm.nb", se=TRUE)+
  labs(y="Invert. Richness", x= expression(paste("Elipsoid Volume cm"^{3})))

plot_nb4 <- ggplot(cafi_coral, aes(x=volume_pg, y=cafi_richness))+
  geom_point()+
  theme_classic()+
  geom_smooth(method="glm.nb", se=TRUE)+
  labs(y="Invert. Richness", x= expression(paste("Software Volume cm"^{3})))

#Create fig. 5 and annotate 
fig5 <- ggarrange(plot_nb1, plot_nb2, plot_nb3, plot_nb4,
          ncol = 2, nrow = 2)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

annotate_figure(fig5,
                top = text_grob("Negative binomial", color = "black", face = "bold", size = 14),
                fig.lab = "Figure 5", fig.lab.face = "bold")

ggsave("fig5.pdf", plot = fig5, path="figs", width=5, height=4)

#Create data tables from negative binomial results for abundance
pl_nb <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid Volume")
pl2_nb <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software Volume")

tab_model(glm.nb1,glm.nb2,
          pred.labels = c(pl_nb,pl2_nb),
          show.se = TRUE,
          show.aicc = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Negative Binomial Abundance"))

	Negative Binomial Abundance
Predictors	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value
Intercept	12.78321	0.15554	9.12736 – 18.08719	<0.001	12.05323	0.15249	8.85505 – 16.55706	<0.001
Elipsoid Volume	1.00004	0.00001	1.00001 – 1.00007	0.002
Software Volume					1.00009	0.00003	1.00003 – 1.00015	0.001
Observations	29				29
R2 Nagelkerke	0.283				0.414
AIC	218.898				215.597

#Create data tables from negative binomial results for richness
pl3_nb <- c(
  '(Intercept)' = "Intercept",
  volume_field = "Elipsoid Volume")
pl4_nb <- c(
  '(Intercept)' = "Intercept",
  volume_pg = "Software Volume")

tab_model(glm.nb3,glm.nb4,
          pred.labels = c(pl3_nb,pl4_nb),
          show.se = TRUE,
          show.aicc = TRUE,
          string.ci = "Conf. Int (95%)",
          string.p = "P-Value",
          digits=5,
          dv.labels = c("Negative Binomial Richness"))

	Negative Binomial Richness
Predictors	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value	Incidence Rate Ratios	std. Error	Conf. Int (95%)	P-Value
Intercept	6.39297	0.10897	5.13025 – 7.92216	<0.001	6.10986	0.10902	4.91180 – 7.54485	<0.001
Elipsoid Volume	1.00002	0.00001	1.00001 – 1.00004	0.004
Software Volume					1.00006	0.00002	1.00002 – 1.00009	0.001
Observations	29				29
R2 Nagelkerke	0.313				0.412
AIC	155.154				152.683

#Probability that the "ith" model minimizes info loss
library(AICcmodavg)

## Warning: no function found corresponding to methods exports from 'raster' for:
## 'wkt'

aicc_field <- AICc(glm.nb1) #calculate AICc from field and software models
aicc_pg <- AICc(glm.nb2)
#AIC_min
p=exp((aicc_pg-aicc_field)/2)
p

## [1] 0.1919075

With \(p=0.108\) and \(p=0.233\), we fail to reject the null hypothesis that the two models comparing elipsoid and software volume are indistinguishable for both abundance and richness using the Vuong test comparing non-nested models.

Theoretically, the lower AIC values for software estimated volume estimating log-likelihood as shown in the tables reflect a relatively better quality model to fit the data. Because a likelihood ratio test requires a nested model assumption, we used the equation \(p=exp^{AIC_{min} - AIC_i/2}\) to calculate the relative likelihood of the two models. Based on our calculation, the invertebrate model with elipsoid volume as the predictor variable was 0.19 times as probable as the model with software volume to minimize information loss.